Quantum fluctuations around black hole horizons in Bose-Einstein condensates 
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We study several realistic configurations making it possible to realize an acoustic horizon in the 
flow of a one dimensional Bose-Einstein condensate. In each case we give an analytical description of 
the flow pattern, the spectrum of Hawking radiation and, the associated quantum fluctuations. Our 
calculations confirm that the non local correlations of the density fluctuations previously studied 
in a simplified model provide a clear signature of Hawking radiation also in realistic configura- 
tions. In addition we explain by direct computation how this non local signal relates to short range 
modifications of the density correlations. 



I. INTRODUCTION 

During the last decade, it has been realized that Bose- 
Einstein condensates (BECs) were promising candidates 
for producing acoustic analogs of gravitational black 
holes, with possible experimental signature of the elusive 
Hawking radiation. The acoustic analogy had been pro- 
posed on a general setting by Unruh in 1981 [l|, and its 
specific implementation using BECs has been first pro- 
posed by Garay et ai, followed by many others [1]. We 
are now reaching a stage where experimental realizations 
and the study of these systems are possible and it is 
important to propose realistic configurations of acoustic 
black holes and possible signatures of Hawking radiation 
in BECs. Other experimental routes for observing analog 
Hawking radiation effects are based on non linear opti- 
cal devices Q or surface waves on moving fiuids Q . Note 
that this last option is restricted to the stimulated regime 
where the Hawking radiation results from a disturbance 
external to the system. 

In this line, density correlations have been proposed in 
Ref. [6] as a tool making it possible to identify the spon- 
taneous Hawking signal and to extract it from thermal 
noise. The physical picture behind this idea is the same 
as the one initially proposed by Hawking 0, Q: quan- 
tum fiuctuations can be viewed as constant emission and 
re-absorption of virtual particles. These particles can 
tunnel out near the event horizon and are then separated 
by the background fiow (which is subsonic outside the 
acoustic black hole and supersonic inside), giving rise to 
correlated currents emitted away from the region of the 
horizon. In contrast to the gravitational case, the experi- 
mentalist is able to extract information from the interior 
of an acoustic black hole. It is thus possible to get in- 
sight on the Hawking effect by measuring a correlation 
signal between the currents emitted inside and outside 
the black hole. This two-body correlation signal appears 
to be poorly affected by the thermal noise and seems to 
be a more efficient measure of the Hawking effect than 
the direct detection of Hawking phonons (see Ref. 9]). 

In one dimensional (ID) flows of BECs, following a 
suggestion by Leonhardt et al. [lo| . it is possible, within 
a Bogoliubov treatment of quantum fiuctuations, to give 



a detailed account of the one and two-body Hawking sig- 
nals. This idea was fully developed in Refs. @ and hj 
to obtain physical predictions for specific configurations. 
Ref. focused on a schematic black hole configuration 
introduced in [6*1 and denoted as a "fiat profile configu- 
ration" in the following: it consists of a uniform fiow of 
a ID BEC in which the two-body interaction is spatially 
modulated in order to locally modify the speed of sound 
in the system - forming a subsonic upstream region and 
a supersonic downstream one - although the velocity and 
the density of the flow remain constant. However, this 
type of fiow, with a position dependent two-body interac- 
tion allowing an easy theoretical treatment, is only possi- 
ble in presence of an external potential specially tailored 
so that the local chemical potential remains constant ev- 
erywhere (see details in Sec. Ill Ap . This makes the whole 
system quite difficult to realize experimentally. 

In the present work we propose simpler sonic analogs 
of black holes for which a fully analytic theoretical treat- 
ment of the quantum correlations is still possible. We 
present a detailed account of the Bogoliubov treatment 
of quantum fiuctuations in these settings and show that 
density correlations provide, also in these realistic con- 
figurations, a good evidence of the Hawking effect. We 
also discuss the recent work of Franchini and Kravtsov 
(1^ who proposed an interesting scenario for explaining 
the peculiarities of the two-body density matrix 

in 

presence of an horizon. Elaborating on the similarities of 
g^^-* with the level correlation function of non standard 
ensembles of random matrices [isj . one can argue that 
the non local features of g^"^^ typical for Hawking radi- 
ation should be connected to a modification of its short 
range behavior. We spend some time for precisely dis- 
cussing this point in the framework of the Bogoliubov 
description of the fiuctuations. Our analytical study of 
the wave functions of the excitations makes it possible to 
obtain an non-ambiguous confirmation of this hypothe- 
sis. 

The paper is organized as follows. In Sec. |n]we present 
three configurations allowing to realize an acoustic hori- 
zon in a ID BEC. Then, in Sec. IIIIl we discuss the practi- 
cal implementation of the Bogoliubov approach to these 
non uniform systems. It appears convenient to describe 
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the behavior of the excitations in the system in terms of a 
^-matrix, the properties of which are discussed in detail. 
This allows to describe the system using an approach 
valid for all possible black hole configurations. Within 
this framework, we study in Sec. IIVI the energy current 
associated to the Hawking effect and in Sec. |V]the den- 
sity fluctuations pattern, putting special emphasis on its 
non local aspects. As discussed above we consider in de- 
tail their connection to short range modifications of the 
correlations. Finally we present our conclusions in Sec. 
IVII Some technical points are given in the appendices. In 
Appendix |X] we present the low energy behavior of the 
components of the S'-matrix, in Appendix [B] we derive an 
expression for the energy current associated to the Hawk- 
ing radiation and in Appendix [C] we precisely check that 
the two-body density matrix fulfills a sum rule connect- 
ing the short and long range behavior of the correlations 
in the system. 

II. THE DIFFERENT BLACK HOLE 
CONFIGURATIONS 

We work in a regime which has been denoted as "ID 
mean field" in Ref. 14]. In this regime the system 
is described by a ID Heisenberg field operator ^{x,t), 
solution of the Gross-Pitaevskii field equation. Writing 
^{x,t) — <t{x,t) exp{—ifit/h) this reads 

ihdt^^~ — dl'l> + [U(x)+gn~^j]'S>. (1) 
2m 

In Eq. ([I} /X is the chemical potential, fixed by boundary 
conditions at infinity; h{x,t) — <I>t<I> is the density oper- 
ator and U{x) is an external potential (its precise form 
depends on the black hole configuration considered), g is 
a non linear parameter which depends on the two-body 
interaction within the BEG and on the transverse con- 
finement. Both are possibly position dependent. For a re- 
pulsive effective two-body interaction described by a pos- 
itive 3D s-wave scattering length a and for a transverse 
harmonic trapping of pulsation uj± , one hasg = 2ahujj_ 
[l5l | . In the flat profile configuration of Ref. p g depends 
on the position x (see Sec. Ill A|) . whereas it is constant 
in the realistic configurations introduced below and re- 
spectively denoted as delta peak (Sec. Ill Bp and waterfall 
(Sec. Ill Cp configurations. 

Within the Bogoliubov approach, in the quasi- 
condensate regime, the quantum field operator $ is se- 
parated in a classical contribution $ describing the back- 
ground flow pattern plus a small quantum correction ip. 
In all the configurations we consider, the flow pattern is 
stationary and one thus writes 

^x,t) = <^{x)+i!{x,t), (2) 

<i?(x) being the solution of the classical stationary Gross- 
Pitaevskii equation 

li'f = -^dl^+[Uix)+gm^]^. (3) 



A black hole configuration corresponds to a disymme- 
try between the upstream flow and the downstream one, 
separated by the event horizon. In the following we use 
a subscript "u" for upstream and "d" for downstream. 
The downstream region corresponds to a: > and is su- 
personic. The upstream region corresponds to a; < and 
is subsonic (see, however, the remark at the end of Sec. 
IIIGp . We thus write 

\ \Aidexp(ifcdx)(/)ti(x) for a; > 0. 

In Q lim^j^-oo |'/'ii(a;)| = 1 and lima;^+oo \4>d{x)\ = 1, so 
that n„ and rid are respectively the upstream and down- 
stream asymptotic densities. Also = mVa/h {a — u 
or d) , where Vu is the asymptotic upstream flow velocity 
and Vd the asymptotic downstream one (Vu and Vd are 
both positive). 

In the following we denote the asymptotic velocities 
of sound as c„ and Cd with mc^ = gana, where gu,d = 
limj;^_oo^+oo g{x) (we keep the possibility of a position 
dependent g coefficient in order to treat the flat pro- 
file configuration of Ref. @). We also introduce the 
healing lengths = h/(mca) and the Mach numbers 
Mo, — Va/ca- In a black hole configurations m^, < 1 and 
md > 1. 

Denoting Uu,d — lim2;_j._oo,+oo U (x) one gets from ([3]) 
and (HI 

-z-^ + Ua+ gana = M and n„K = naVd- (5) 
2m 

The first of these equations corresponds to the equality 
of the asymptotic chemical potentials and is required for 
a stationary flow; the second equation corresponds to 
current conservation in a stationary Sow. 

The precise form of the flow pattern is specifled by the 
functions (puix) and 4>d{x) which depend on the configu- 
ration considered. In all the configurations treated below 
(j)d{x) is a constant of the form 

(l)d{x) ^ cxp(i/3<j), (6) 

meaning that the downstream flow pattern is flat with 
a constant density and velocity. The value of fid de- 
pends on the configuration considered. As for the up- 
stream flow pattern, the stationary flow condition im- 
poses \\mx-^-oo <Puix) — exp(i/3„), where /Su is a constant. 

After having defined the notations and the common 
aspects of all the flow patterns, we now give the precise 
value of the configuration-dependent parameters. 

A. Flat profile configuration 

We recall here the value of the parameters in the flat 
profile configuration studied in 0, Q . In this case the (j)a 
functions of Eq. (|4]) assume a very simple value: (j^uix) = 
4>d{x) — 1 (and thus /3„ = /3d = 0). One has 

U(x) = <f < °' (7) 

\ Ud for X > 0, ^'^ 
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and 



9u 
9d 



for a; < 0, 
for x > 0, 



(8) 



chosen so that a flow with Vu — Vd = Vq and = = 
no is solution of Eqs. (|3]) and ([5]); i.e., Eq. (|4]) reduces 
to = y^TiQ exp(ifco2^) for aU x (fco = mVo/h). This 

imposes 



and 



g„no + Uu= gdriQ + J7rf. 



(9) 



(10) 



We finahy note that in the flat profile configuration one 
has Cd<Vd^Vu < Cu- 

In the numerical simulations of Rcfs. H, [ll|, a genera- 
lization of this step-hke configuration has been used; one 
considers smooth U{x) and g{x) functions imposing the 
continuous version of ([TU| : g(x)nQ + U{x) = C^'. The 
theoretical approach is the same as in Ref. [9] but the 
Bogoliubov-de Gennes equations [Eq. (IT9l) below] have 
to be solved numerically whereas the step-like configu- 
ration characterized by Eqs. ([7|) and © allows for an 
analytical treatment. 

The flat profile configuration can be numerically im- 
plemented in a dynamical way as explained in Ref. Q. 
However, it is fair to say that the corresponding experi- 
ment seems rather difficult to realize. Moreover the flat 
profile configuration is very sensitive to the total atom 
number, a quantity which is not easily controlled experi- 
mentally. Besides a local monitoring of g{x) has not yet 
been demonstrated. There are the reasons why in the fol- 
lowing subsections we introduce two new types of sonic 
horizon which can be implemented experimentally more 
easily. 



B. Delta peak configuration 

In this configuration the non linear coefficient g is con- 
stant and the external potential is a repulsive delta peak: 
U(x) = A6{x), with A > 0. It has been noticed in Ref. 
[16| that one can find in this case a stationary profile 
with a flow which is subsonic far upstream and super- 
sonic downstream (i.e., a black hole conflguration). The 
upstream flow corresponds to a portion of a dark soliton 
profile. More precisely, for x < 0, one has 



tanh 



Xo 



cos 9 ] — i sin ( 



(11) 



where sinfl = m„, and one can restrict oneself to 9 € 
[0, tt/2] (then = tt+9). As is also the case for the other 
configurations studied in the present work, the down- 
stream flow has a constant density and velocity [cf. Eq. 
([5])]. The typical profile is displayed in Fig. [T] 



U{x) = AS{x] 
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n{x) 



Vd > Cd 
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FIG. 1. (Color online) Density profile in the delta peak con- 
figuration. The flow is directed toward positive x. The delta 
potential is represented by a (red) vertical straight line. The 
density in the upstream region {x < 0) is a portion of a dark 
soliton (see the text). The region a; > is supersonic. It 
is shaded in the plot for recalling that it corresponds to the 
interior of the equivalent black hole. We keep this convention 
in Figs. [1 [3] and H 



Once m„ = Vu/cu is fixed (< 1) all the other parame- 
ters of the flow are determined by Eqs. ([S|). Defining 

2 I one gets 



y 



(-l + v/l + 8/«.2) 



rid 



Vd 

Vu 



md 
mu 



Cd 
Cu 



1 



(12) 



By imposing continuity of the wave function [$(0) — 
Y^rid exp(i/3(;) = ^/nZ 't'uiO)] and the appropriate mat- 
ching of its first derivative [dx^iO~^) — dx^iO~) = 
2mft~^A<i>(0)] one gets 



sin /3d = -muy/y, 



Xq 



cos 9 



tanh 



ly-1 



tan^? 



and also 



A 



with A 



y-1 



(13) 
(14) 

(15) 



In this configuration one has Vu < Cd < Cu < Vd which 
corresponds to a black hole type of horizon. Related work 
for a double barrier configuration recently appeared in 

Note that the flow depicted in Fig. [1] corresponds to 
a very speciflc case in the parameter space spanned by 
the intensity of the delta potential and the flow veloci- 
ties, which is on the verge of becoming time dependent 
(see Ref. fl^). This is reflected by the fact that, for this 
configuration, m„ and md cannot be fixed independently 
(in contrast to what occurs for the flat profile case) . One 
might thus legitimately expect to face a fine tuning prob- 
lem to experimentally fulfill all the required boundary 
conditions ([H]) and 1^^. One could also ar- 

gue that a delta potential is non standard and that the 
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FIG. 2. (Color online) Same as Fig. [T]for the waterfall con- 
figuration. 



specific structure of the flow pattern displayed in Fig. [T] 
would disappear for a more realistic potential. However, 
it is shown in Ref. that this configuration can be 
rather easily obtained by launching a ID condensate on a 
localized obstacle (not necessarily a delta peak). In this 
case, there exists a sizable range of parameters where, 
after ejection of an upstream dispersive shock wave, the 
long time fiow pattern becomes of the type illustrated in 

Fig.m 



C. Waterfall configuration 

In this configuration the two-body interaction is cons- 
tant and the external potential is a step function of the 
form U (x) — —Uo'd{x), where 8 is the Heaviside function 
(and Uq > 0). In this case, a stationary profile with a fiow 
which is subsonic upstream and supersonic downstream, 
i.e., a black hole configuration, has been identified in Ref. 
[igj . The upstream profile is, as for the delta peak con- 
figuration, of the form (jlip . with here xq = 0, i.e., the 
upstream profile is exactly one half of a dark soliton. The 
corresponding density profile is displayed in Fig. [31 

The equalities ^ and the continuity of the order pa- 
rameter at the origin impose here 



exp(i/3d) 



rid 



1 



-i, and 



gnu 



(16) 



(17) 



In this configuration one has Vu ~ Cd < Cu < Vd which 
corresponds to a black hole type of horizon. 

The remark given at the end of Sec. Ill Bl is here also 
in order: there might be a fine tuning problem for ver- 
ifying Eqs. din]) and pT| . Although we did not per- 
form here the time dependent analysis done in Ref. [Tsj 
for the delta peak configuration, we believe that, also 



in the present case, it is possible to dynamically reach 
the stationary configuration depicted in Fig. [51 This is 
supported by the experimental results presented by the 
Technion group [s*] who studied a very similar configura- 
tion (with the additional complication of the occurrence 
of a white hole horizon). These results show no impor- 
tant time dependent features near the black hole horizon 
and we are thus led to consider that the stationary wa- 
terfall configuration of the type illustrated in Fig. [2l is 
stable and can be reached experimentally. 

We make here a remark which is also relevant for the 
delta peak configuration: the precise location of the sonic 
horizon is not well defined. In both configurations (water- 
fall or delta peak) one may define a local sound velocity, 
and the point where the local fiow velocity exceeds the 
local speed of sound can be chosen as the location of the 
sonic horizon. Then one finds that the sonic horizon is 
located slightly upstream the interface x = (in the wa- 
terfall configuration for instance, at a; = the local fiow 
velocity is already times larger than the local sound 
velocity). However, the local sound velocity is an appro- 
ximate concept, only rigorously valid in regimes where 
the BEC density varies over typical length scales much 
larger than the healing length. This is not the case in the 
waterfall and delta peak configurations near x = and, 
as a result, the concept of sonic horizon is ill defined. In 
the fully quantum treatment presented below we do not 
use this concept: the important point for our analysis is 
simply that the upstream fiow velocity is asymptotically 
(i.e., when x — > — oo) larger than c„. Hence, for precise- 
ness we do not state that the upstream flow is subsonic, 
but that it is asymptotically subsonic. 



III. FLUCTUATIONS AROUND THE 
STATIONARY PROFILE 

In this section we establish a basis set in each of the 
flow regions (upstream and downstream) which will be 
used in Sec. IIII FI for describing the quantum fluctuations 
in the system. The simplest way to obtain this basis set is 
to start from expression with $(a;) given by ([l]), and 
to treat tl){x,t) as a small time dependent classical field, 
denoted as il){x,t) in the present section, with <^{x) + 
ip{x,t) solution of the classical version of ([T]). One looks 
for a normal mode of the form 



V'(x,i) 



AkaX 



(x, uj)e' 



-iujt 



- w, 



,(x,a;) 



(18) 



with a = u for x < and a = d for a; > 0. $(a;, t) defined 
by Eqs. ([2]) and (fT8)) describes small oscillations with 
pulsation oj of the order parameter around the ground 
state ^{x). In the following we drop the oj dependence of 
functions Ua and Wa for legibility. We also write Xa = 
xj^a (and then k^x — maXa)- Linearizing the Gross- 
Pitaevskii equation, one gets at first order in 



— C 



(19) 
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with 

where = f^/{gana) and Ha = + 2\(j)a\'^ — 

1. Hence, the column vector formed by Ua and Wa is 
an eigen-vector of the so caUed Bogohubov-de Gennes 
Hamihonian Ca- 

The present section is organized as foUows. We first 
consider solutions of (US]) for e M in Sec. IIII Al We 
give the expression of these solutions in sections [ill Bl and 
IIII C[ specifying only what we need for the following step: 
that is, for a = d, we only display the form of the solution 
when a; > 0, and for a = u, we only display the form of 
the solution when a; < 0. The most general fluctuation of 
pulsation is a linear combination of eigen- modes for the 
upstream region glued at a; = with a linear combination 
of the downstream eigen-modes. We explain how this 
matching is done in Sec. IIII Dl Finally, in Sec. IIII El we 
specify the form of the scattering modes which are the 
appropriate modes used for quantizing the fluctuations 
in Sec. IIITfI 

A. Properties of the eigen-functions of the 
Bogoliubov-de Gennes equation 

The relevant eigen-functions of (IT^ are of the form 

[Mx))~ [mix))' ^^'> 

where the functions Ui(x) and YViix) are constant for 
|a;| — ?> CXI (more precisely in the domain where is con- 
stant). Their exact form will be specified later [Eqs. ([?7)) 
and ([32]) ]. In Eq. ((2T|) the Qis are the dimensionless 
wave vectors of the Bogoliubov modes, solutions of 

{Sa-maQf ^LoliQ), (22) 

where 



WB(Q) = Qyi + ^ (23) 

is the Bogoliubov dispersion relation in a condensate at 
rest (written in dimensionless form). Note that Qi ~ 
solution of (|22]) - is sometimes complex; this fact is taken 
into account in the following. In particular, for a — u 
{a = d) one should discard values of the wave vector 
such that lai{Qg) > {lm{Q£) < 0). This corresponds to 
eliminating the evanescent channels in the region where 
they are divergent. For instance, a mode with lm(Qi) < 
diverges when x — >■ -|-oo, which will correspond to the 
supersonic region (labeled d) in the following, and we 
thus discard it. The dispersion relations and the different 
real wave vectors are displayed in Fig. [3] 

The index a of Eq. ^51) is specified to £ in (pij) for 
identifying the branch of the dispersion relation to which 
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FIG. 3. (Color online) Dispersion relation (|22p . In each plot 
the horizontal dashed line is fixed by the chosen value of ui. 
The gf(a;)'s are the corresponding abscissae. Only the real 
eigen-modes are represented; their denomination is explained 
in the text; their direction of propagation (left or right) is 
represented by an arrow. The part of the dispersion relation 
corresponding to negative norm states (see the text) is repre- 
sented with a dashed line. The upper plot corresponds to a 
subsonic flow. The lower one corresponds to a supersonic flow; 
it is shaded in order to recall that it describes the situation 
inside the black hole. 



the considered excitation pertains. For precise notations, 
t is taken as a double index, because it is clear from (|22|) 
that the values of the wave vectors are not the same in 
the subsonic and supersonic regions, i.e., they depend on 
a. Specifically, when a — u, I & {u|in, u|out, u|eva}. 
When a — d^ there are two cases, depending if uj is 
lower or greater than a certain threshold VL. If a; < il, 
i e {dl|in, dl|out, d2|in, (i2|out} and when uj > VL, I e 
{dljin, dljout, djeva}. 

We have chosen to label the real eigen-modes as "in" 
(such as dljin for instance) or "out" (such as u|out) de- 
pending if their group velocity [its explicit expression is 
given below, Eq. (j29p ] points toward the horizon (for 
the "in" modes) or away from the horizon (for the "out" 
modes) in a black hole configuration, i.e., with a sub- 
sonic region at left of the horizon, the supersonic region 
being at the right. The wave vectors labeled wjeva and 
d|eva are complex and correspond to evanescent channels 
(as explained above, one selects the complex Qf's which 
describe waves decaying at infinity). 
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The threshold appearing in the lower plot of Fig. [3] 
is reached only for a supersonic flow, for a wave vector 
q* such that 



+ f + Y 



(24) 



The existence of this threshold is a consequence of the 
behavior of the large momentum part of the dispersion 
relation uJb{Q) of a condensate at rest. More precisely, 
the part of the dispersion relation which shows a local 
maximum in the supersonic region corresponds to the 
particular solution of (1221) where — m^^Q = —uj^{Q). 
At Q — Q*^ one has exactly dajg/dQ = m^, and for 
Q > Q*i, one has dus/dQ > md- Hence the part of 
the spectrum with Q > Q*^ corresponds to excitations 
whose group velocity in the frame where the condensate 
is at rest {dajg/dQ) is larger than the flow velocity md 
(we use here dimensionless quantities). One can have 
duJe/dQ > md > I only because the dispersion relation 
(1^5)1 grows faster than linear at large Q. We will see in 
Sec. IIVI that the corresponding waves play an important 
role in the zero temperature Hawking signal. 

It is easy to verify that, if (uq, iDq) is a solution of 
associated to an eigen-energy e^, then (w* , u* ) is also a 
solution of ([T9l). now associated to eigen-energy — £«. Be- 
sides, from expression (jl8p . one sees that both solutions 
describe the same perturbation of the condensate. As a 
result, one can always select eigen-modes with positive 
values of = hu)/{gana), and in all the following we 
chose cj € 

One may also notice that the above symmetry of the 
wave function is associated to the normalization of the 
eigen-modes. For instance, it is clear that the symme- 
try operation changes the sign of jW^p — jW^P and one 
can show by simple algebraic manipulations that, for real 
Qe, the sign of — jW^P is the same as the sign of 
£a — maQe- In the following we denote the eigen-modes 
for which \Ue\'^ — |W£p > as having a positive norma- 
lization (for a recent discussion of this point see see, e.g. 
Ref. asweUasRefs. [21,22]). In Fig. [3] their disper- 
sion relation is represented with a solid line, whereas the 
eigen-modes with negative normalization are represented 
with a dashed line. In particular, in the supersonic re- 
gion, the c?2|in and d2|out channels have negative norm 
for < w < 

It was shown in [2^ that each eigen- vector of equation 
([TO)) is associated with a conserved (i.e., x independent) 
current Jg. We show below (see Sec. IIV[) how Jg relates to 
the energy current in the system. Ji is zero for complex 
Qi (evanescent waves do not carry any current); this can 
be proven directly in our specific case, but we do not 
display the proof here. For real Qe one gets 

Jl = Ca {Qe + ma)\Uf]'^ + {Qi ~ mo,)\yVe\^ 

+ c„Im {Uldx^Ue + W^dx^ Wi) . (25) 
Going back to dimensioned quantities and using the ue 



and We functions, this reads 
h r 



2m 



u*e{ka-idx)ui-wg{ka + idx)we +c.c., (26) 



where "c.c." stands for "complex conjugate" . 



B. Downstream region: a; > 

We recall that in this region the flow is supersonic 
and that the eigen-vectors are labeled with an index 
i e {dl|in, dl|out,d2|in,d2|out} when uj < D, and £ G 
{(il|in, dljout, djeva} when uj > fl. 

Here (f)'^ appearing in Eq. (1^ does not depend on 
X and is equal to exp{2i/3d) 124!]. This implies that the 
functions Ue and We are also x independent. One finds 



Ue 

m 



Ce 



- Eeyl^'' 
Ee)e-'^'' 



(27) 



with Ee — £d — mdQe and Ce is a normalization constant 
which we always chose real and positive. The correspond- 
ing current is easily evaluated using Eq. ([25]) . For real 
Qe one gets 



where 



Ji = Vg{Qe){\Ue\^ - jW^H, 



1/ fn \ 



(28) 



(29) 



is the group velocity in the laboratory frame and qe — 
Qi/Ca with here a = d, but Eqs. (1^ and (1^ are valid 
also for a ~ u. 

A typical choice for the normalization constant is Ci = 
|2Re(i;;Q2)|i/2, This ensures that \Ue\^ - \yVe\'^ = ±1. 
For our case, it is more appropriate to multiply the pre- 
vious expression of Ce by |V^(Q£)|"'^^^, so that for real Qe 
one has Je = ±1. Hence we chose 



Ce = |2Re(^; Q2)^^(g^ 



which implies 



±1 



\Vg{Qi)\' 



(30) 



(31) 



The modes for which the factor +1 appears in the above 
expression are the positive norm modes previously dis- 
cussed. The others are the negative norm modes. With 
the normalization (pij) . from Eq. (pS)) . one sees that 
Je = 1 either for non evanescent modes of positive norm 
propagating to the right, or for non evanescent modes of 
negative norm propagating to the left. In the other non 
evanescent cases (modes of positive norm propagating to 
the left or modes of negative norm propagating to the 
right) Je = -1. More precisely: Jdi|out = +1 = Jd2\in 



and Jdi\in = —1 — Jd2\out- We will see below (Sec. |IV 
that hujJe is the energy current associated to mode £ and 
thus negative norm modes can be interpreted as carrying 
negative energy. 
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C. Upstream region: a; < 

We recall that in this region the flow is asymptot- 
ically subsonic and that i G {wjin, u|out, w|eva}. In 
the flat profile configuration one has = 1 and 

the functions Ui and have the same form as the 
ones displayed in the previous section. Hence in the 
remainder of the present subsection we concentrate on 
the delta peak and waterfall configurations where de- 
pends on Xy^. In this case the functions Ui and have 
a more complicated expression than in the downstream 
region (see, e.g., Appendix A of Ref. [25|)- Defining 
x(X„) = cos0tanh[(A„ — Xq) cos6')], where ATq = xq/^u 
(we recall that in the waterfall configuration = 0), one 
gets 



1 



[Q«/2 + e„/Q, + ix(A„)]' 
[Q,/2 - e„/Q, + ix(X«)]' 



(32) 



where TD^ is an arbitrary constant, the value of which is 
determined by the normalization (see below). The cur- 
rent (^5]) corresponding to Ut and given in ([5^ is 
most easily evaluated at — )■ — cxi, i.e., in a region 
where Ui and become independent of Xu- In this 
region one has 



1 



[Qi/2 - e^/Qi 



ICOSf 

icos( 



(33) 



Since tends to a constant [exp(i/3„)] when X^ — > — oo, 
one could, in this region, use for the Bogoliubov modes 
an expression similar to Eq. (j27|) which is used in the 
downstream domain. Indeed, it is difficult to see it from 
the above formula, but we have checked that (1551) is pro- 
portional to an expression similar to (|27p where Pa is 
replaced by Pu- However, expression p3|) is here more 
appropriate since its position-dependent version ([5^ is 
valid for all X„ < 0. From ([55)) one gets for real Qg 



\uif - \m 



V Q ^ 

>-oo 12?, I 



Qe 



(34) 



where Ef = Eu — muQi. In the following, the constant 
will be chosen such that = ±1 for real (see 
the discussion at the end of Sec. IIII B| one has here 
Ju\out = —1 and Ju\in = +!)■ Also, we chose I?„|in on the 
positive imaginary axis in the complex plane and | out 
on the negative imaginary axis. This implies that, for 
real Qi, 



\Qe\ 



(35) 



The particular choice of phase in ([55]) is based on aes- 
thetic grounds: it ensures that the w — >■ limit of the 
eigen-function (j55]) upstream the horizon in the water- 
fall and delta peak configurations has the same phase as 
its equivalent for the flat profile configuration. This will 



make it possible in Appendix [Cl to obtain formulae valid 
for all three types of configurations [Eqs. (|C4p and dCSp ]. 

If Qe is complex, the expression is more complicated; 
we write it here for completeness. One takes 

Ve = \VgiQe)\^ x 8Re[Ei{eu/Qi)^] 



4 Eu cos^ 9 



2i £„ cos 6 



\Qi\' 

Qi — Q*e 



(36) 



This expression is clearly real and it ensures that, as in 
the downstream region, the normalization pip is fulfilled 
for all Qi. 



D. Matching at a; = 



Let us denote 



5q(x) 



_ f ex.p{imaXa)Ua{x) \ _ f Ua{x) 



eyi-p{ — lmaXa)Wa{x] 



Wa{x) 



and 



exp(i((5f 
exp(i((5^ 



- ina)Xa)Ui{x) 
m^)X^)Wl{x) 



(37) 



(38) 



Remember that the index a is equal to either u or 
d, depending which side of the horizon one consi- 
ders, whereas (. labels the eigen-modes of Eq. (jl9p : 
I £ {wjin, wjout, wjeva, dljin, dljout, (i2|in, c?2|out, d|eva}. 
More precisely, S„ describes the excitations in the sub- 
sonic region; it is a linear combination of ^u\im ■^ulout 



and ii„ 



which describes the same excitation in 



the supersonic region, is a linear combination of .^dil 



^cZl|out: ' — 'd2\ini "d2|out 



and S, 



Then the matching conditions at the horizon read 



S,(0) = Sd(0), 



and 



2m 



dSd 
dx 



(0) 



dS„ 
dx 



(0) 



= AS„(0). 



(39) 



(40) 



In the case of the fiat profile or of the waterfall configu- 
ration, Eq. (Uni) also holds, but then A = 0. 



E. The scattering modes 



Amongst all the possible modes described in Sec. IIIIDI 
as linear combinations of the S^'s, we are primarily inter- 
ested in the scattering modes. These are the three modes 
which are impinging on the horizon along one of the three 
possible ingoing channels: u|in, dljin or (i2|in. Each of 
these ingoing waves gives rise to transmitted and reflected 



waves which, together with the initial ingoing compo- 
nent, form what we denote as a "scattering mode". It is 
natural to label these modes according to their incoming 
channels, but since each mode includes more than the 
ingoing wave that generates it, for avoiding confusion in 
the notations, we use capital letters and denote the scat- 
tering modes as S", 5°"^ and S"^. 

For concreteness, we now give the expression of each 
of the scattering modes. Each mode has a different ana- 
lytical expression on each side of the horizon. According 
to our conventions we denote these expressions as S^(a;), 
S°^(a;) and E°^{x) in the upstream region and SJ^(a:), 
SJ^(a;) and SJ^(x) in the downstream one. Specifically, 
one has 

— Sdl,u^dl\out + — ^)Sd2,u'^d2\out 

'^Di _ Q -ji^ I qcva -izr 



U mode (initiated by u\m) 



(41) 



^d^ = '^dllin + 'S'dl.dl'^dllout 

+ 8(r2 — Uj)Sd2,dl'^d2\ont 

+ q{ll! - ri)S';^y|;^Sd|eva, 

~ ©(^ ^ ^){Su,d2^u\out + S'u^d2^u\cva.) ^ 

= ^ '^)(2d2|in + Sdl,d2'^dl\ovLt 

+ 'S'd2,d2"d2|out)- 



For legibility, the x dependence of the S vectors has 
not been displayed in the equations. The three differ- 
ent modes are displayed in a pictorial way in Fig. 0] 
where the purple wiggly lines correspond to the evanes- 
cent modes wjeva and djeva that cannot be represented 
in Fig. 131 Note that when w > fi, the outgoing d2 wave 
(involved in the U and Dl modes) becomes evanescent 
and the D2 mode disappears, because the incident seed 
for this mode (the propagating d2|in wave) disappears. 
This is taken care of in formulae (PT|) by the Heaviside 
functions 0(r2 — w) and <d{uj — U). 

The S'-coefficients (S'„,„, S'°y^, etc.) in Eqs. (gl]) are 
complex and do not depend on x (they do depend on lo 
though). For each of the three scattering modes, one has 
four such coefficients which, once the incident channel is 
fixed, are determined by solving the 4x4 system of linear 
equations ([M)) and pOj) : hence, the 5- parameters depend 
on the configuration considered (flat profile, delta peak 
or waterfall). Physically, the square moduli |5'^'_^(tLi)p 
of the S'-matrix elements give the transmission or reflec- 
tion coefficients for a i/-ingoing mode of energy fvx) which 
scatters into an z^'-outgoing mode at the same energy. 

Current conservation can be written in a simple ma- 
trix form provided the normalization of the real modes 
is defined in such a way that J( = ±1 (as done in Sees. 
IIII Bl and IIII Cp . Defining, for w < 57, the 5-matrix as 



u 


u|in 
out w|eva 


dl out d2 out 
(d eva) 


Dl mode (initiated by dl in) 


u 


out w|eva 


dl out d2 out 
dl in (d eva) 


D2 mode (initiated by (i2 in) 


u 


out u 1 eva 


dl out d2 out 
^^d2|ir^* 



FIG. 4. (Color online) The scattering modes. The color code 
is the same as in Fig. [S] The additional purple wiggly lines 
correspond to evanescent channels. For uj > i}, the D2 mode 
disappears and the channel d2|out is replaced by d|eva in the 
two upper panels {U and Dl modes). The region correspond- 
ing to the interior of the black hole is shaded as in the previous 
figures. 



current conservation reads 

S^T]S = Tj = St]S\ where 77 = diag(l, 1, -1). (43) 

The coefhcients such as are not involved in current 
conservation since the evanescent waves carry no current. 
Note that for ui > fl, the iS-matrix is 2 x 2 because, 
the outgoing d2 mode - being evanescent in this case - 
is not involved in current conservation: one simply has 
the usual unitarity relation S^S = diag(l,l). We have 
checked that our results for the scattering matrix indeed 
fulfill the 77-unitarity condition for a; < (and the 
unitarity condition for uj > Q). 

In the following we will need to determine the low- 
u behavior of the components of the ^-matrix. In the 
three configurations we considered, we always find that, 
for v ~ u.dl, d2, one has 



S,,di ^^ + h,MiV^+0{el/^), 
S,M2 ^^ + K,d2V^ + 0{el/^), 



Su,u Su,dl Su,d2 
S = I Sdl,u Sdl,dl SdlA2 
^Sd2,u Sd2,dl Sd2,d2^ 



(42) 



(44) 



where e„ = huj/{mc^) and the /'s and the /I's are di- 
mensionless complex numbers. We determined them an- 
alytically in the three configurations we considered. The 
relevant formulae are given in Appendix [Xj 
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F. Quantization 

The field operator ip{x, t) associated in the Heisenberg 
representation to the elementary excitations on top of the 
background [as defined by Eq. ([5])] is expanded over the 
scattering modes: 



-'o V27r L 



X, uj e 



*aL(w) 



AljJ 



+ <(x,a;)e'"*aiM 
UD2(a;,a;)e"'"*a|32(a;) 
+ i(};,(x,w)e'"*a„,Hl, (45) 



where we have written explicitly the uj dependence. The 
dt(a;)'s create an excitation of energy hu> in one of the 
three scattering modes ([/, D\ or D2). They obey the 
following commutation relation: 



(46) 



From expression (|^5|) one sees that the D2 mode (which 
originates from the negative norm d2\m. channel) is quan- 
tized in a non standard way: the role of the creation and 
annihilation operators is exchanged compared to the U 
and Dl modes. Using the current conservation relation 
(|43p . one can show that this choice of quantization is 
necessary for fulfilling the appropriate Bose commuta- 
tion relation of the if) operator: 



['4}{x,t),'^\x\t)] ^5{x-x'). 



IV. RADIATION SPECTRUM 



(47) 



The Hawking signal corresponds to emission of radi- 
ation from the interior toward the exterior of the black 
hole . In our specific case the energy current associated 
to emission of elementary excitations is (cf. [26]) 



n(x,t) 



2m 



dt^\x,t) dA{x,t) +\v.c., 



(48) 



where "h.c." stands for "hermitian conjugate". From 
expressions ([2]) and ()45|) one can write the average current 
n(x) = (n(a;,i)) under the form 



/■°° dw ^ , 
Jo 27r 



(49) 



where J{x, co) [and accordingly n(x)] can be separated in 
a zero temperature part Jo(x,aj) [Ho (a;)] and a "thermal 
part" Jt(x,uj) [nr(x)] with 

Jo{x,u}) = ^ - ^ wl{ka+idx)wr. 



and 



where 



Jrix, to) 



+ 9(ri - w)M*2(fcQ - idx)uo2 



^ Ji(x,cj)ni(w), 



(50) 



(51) 



Jr. = 



h 

2m 



ul{ka-idx)uL-wl{ka+idx)wL +C.C. (52) 



and nL{uj) = (di(a;)ai,(w)) is the occupation number of 
the mode L. Note that in expression ([5T|) the D2 mode 
contributes only for cj < f2. Comparing the expression 
(|5^ with (1^51) one sees that Jj, is the conserved current 
carried by a scattering mode L; it is x independent for 
the stationary flows we consider in the present work. 

Note for avoiding confusion that what we call a "zero 
temperature term" is the contribution to the Hawking 
signal that exists even when the system is at zero tem- 
perature. It will be described below (Sec. IIVB|) by an 
effective radiation temperature Th (the Hawking tempe- 
rature) which is not the temperature of the BEC. 



A. Energy current in a black hole configuration 

For large and negative x (i.e., deep in the subsonic 
region) we show in Appendix |B] that formulae (|4T|) and 
([50)1 yield the very natural result 



H 



— 



Ztt 



(53) 



The minus sign in this equation indicates that the ener- 
gy current is directed toward — oo, as clearly seen from 
Appendix [BJ If one computes the energy current for a 
point deep in the supersonic region (i.e., for x large and 
positive) one gets the same result as (I53|) . in agreement 
with the conservation of the energy flux in a stationary 
configuration. Note that the zero temperature radiation 
Hq vanishes in absence of black hole, as expected. In 
presence of a black hole, the integral gives a finite re- 
sult, corresponding to a Hawking signal emitted even for 
T — 0. This remark, together with the specific form of 
Eq. ([53l) . shows that one needs two ingredients for hav- 
ing a T = Hawking radiation from a black hole: (i) a 
d2|in mode and (ii) a c?2 o m mode conversion, i.e., a 
non zero Su.d2 coefficient. Remember that condition (i) 
is fulfilled only because the dispersion relation in a su- 
personic flow bends down at high q (see Fig. [31 bottom 
panel) , which is a consequence of the non linear behavior 
of the Bogoliubov dispersion relation ([^31) . As discussed 
in Sec. IIII Al the d2 incoming channel corresponds to 
waves whose group velocity in the frame of the conden- 
sate is larger than Vd- It is thus not surprising that these 
"fast" modes are involved in the Hawking radiation since 
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they are able to overrun the flow of velocity Vd and thus 
to escape the black hole. This is clearly a flaw of the BEC 
analogy of gravitational black holes: point (i) is certainly 
not fulfilled in the gravitational case since the group ve- 
locity of photons is a constant (the speed of light). As a 
result the number of ingoing and outgoing channels are 
not equal for gravitational black holes and one cannot 
get a stationary description of the Hawking effect: one 
has to take into account the dynamics of the formation 
of the horizon. 

From Eqs. ([TO and ([STt . the term Jt{x,uj) can be 
rewritten as 



Jt(x, uj) 



71„(w)(l - \Sum\'^) 



noi{i^)\Su.di\'^ - no2(w)|S'„,d2p 
2 



-|[7i„(a;) - 7iDi(a;)]|S'„,di| 



[n„(a;) + no2(w)]|S'„,d2p|- 



(54) 



In absence of black hole, the Su.d2 term in (|54l) disap- 
pears. As a result, at thermal equilibrium, i.e., when 
71,7(0;) = noi(w) is a thermal Bose occupation number of 
the form 



1 



Cxp(^) -1' 



(55) 



one has Jt{x,uj) ~ [this is most easily seen from the 
last expression of Jt in Eq. ((54|)]. This is a very pleasant 
result demonstrating that at thermal equilibrium there is 
no Hawking radiation at all for any type of configuration 
connecting two asymptotically subsonic regions. 

In the case where an acoustic horizon is present, a finite 
temperature configuration may be reached in the manner 
presented in Ref. Q: one branches the black hole con- 
figuration adiabatically starting from a system initially 
uniformly subsonic at thermal equilibrium. This changes 
the dispersion relation in the supersonic part, but not the 
occupation number of the adiabatically modified modes 
(see the discussion in [^). In this case ((54)) yields a finite 
value for the H^, term (because Su,u, n^i and are 
regular at low uj). 



B. Hawking temperature 

The zero temperature radiation Hq as given by 
(I53p corresponds to an emission spectrum given by 
\Su,d2{'^)\'^ ■ Can this be described by an effective tempe- 
rature, i.e., can \Su,d2{<^)\'^ be approximated by a factor 
of the type F x nT^{Ljj), where Th is the effective tem- 
perature of radiation? We address this question in the 
present subsection. 

One could first argue that the addition of the "gray- 
ness factor" F is an unnecessary complication of the fit 
of |5'u,(;2p by a thermal spectrum. Indeed, in most cases. 



F is found to be close to 1, but this is not generally true 
(see the discussion below and the inset of Fig. [S|) and this 
is the reason why we keep a certain degree of "grayness" 
in the present analysis. 

Obviously, the identification of |S'u,d2p with a Bose 
thermal factor can only be approximate because, whereas 
a term such as Ut-^Ill!) is finite for all w G K^, |5'«,d2p 
abruptly cancels for w > Vt. Nonetheless one can try 
to find the best possible approximation by comparing 
the low-w expansion of |S'u.d2p [from Eqs. (011)] with 
FxnT^(a;) = T[k^T^/{fiuj)-\+0{i^)]. This immediately 
yields 



-4Re(/:.rf2/i^- 



,d2) 



and 



^B^H \fu. 



d2\ 



met 



(56) 



The analytical expressions obtained in Appendix |X] for 
fu,d2 and hu,d2 in the different configurations we con- 
sider yield the following estimates of the Hawking tem- 
perature: 



met 




(flat profile). 



(waterfall) . 



(57) 



We do not display here the formula for the delta peak 
configuration because it is too cumbersome. Instead we 
show in Fig. [5] the corresponding curve relating to 
m„ in the delta peak configuration and compare it with 
the results of the waterfall configuration. One first no- 
tices from the figure that Th — ?► when mu ~^ 1: this 
is expected because in this case the horizon disappears. 
One sees also that Th remains finite in the limit )«„ — > 
in both configurations. However, this limit is singular in 
the sense that it corresponds to a very peculiar flow. For 
instance, in the waterfall configuration, the analysis of 
Sec. Ill CI shows that the flow with m^, = is observed for 
a step with C/q — >■ cx), and has downstream a zero density 
and an infinite velocity: the corresponding fiow pattern 
is most probably unreachable. Moreover, one sees from 
the inset of Fig. [5] (and also from the analytical expres- 
sion given in Appendix |X| that in this case F — >■ and 
the expected signal disappears. In the remainder of this 
work we rather consider a typical setting with m^, — 0.5. 
For this value of m„ one gets md — 1.83 and F ~ 0.977 
in the delta peak configuration, md — 4: and F ~ 0.980 in 
the waterfall configuration. 

Once Th is determined through the above low energy 
analysis, one should check, as done, e.g., in Ref. (llj . if 
the approximation of |S'u,d2|^ with a thermal spectrum is 
accurate in the whole emission window a; S [0, SI]. This 
is done in Fig. [5] in the case of the delta peak configura- 
tion. One sees from the figure that the overall agreement 
is quite good. The same good agreement is obtained for 
the waterfall and fiat profile configurations, and this le- 
gitimates the definition of a Hawking temperature in the 
three configurations considered in the present work. This 
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FIG. 5. Normalized Hawking temperature A:bT'h/(w.c^) as a 
function of the upstream Mach number = 14 /cu for the 
delta peak and waterfall configurations. The inset displays 
the grayness factor F as a function of mu for these two confi- 
gurations. 



X 
h 

e 

X 

i-i 

PI 



Co 




0.1 0.2 0.3 

u) in units of Cu/^u 



FIG. 6. (Color online) Radiation spectrum in the delta peak 
configuration. Red curve: 15*11,^2 1^ as a function of the dimen- 
sionless quantity to^u/cu- Black curve: F x nT„(i^). The plot 
is drawn for m„ = 0.5; in this case ksT^/{mcl) ~ 0.128 and 
F ~ 0.977. The difference between |S„,d2p and F x riTg is 
maximum when uj = Q, and is close to 0.06 in this case (for 
the chosen value of m„ one has f2^u/cu ~ 0.369). 



was not a priori obvious because the concept of Hawking 
temperature is of semi-classical origin (see, e.g., [13]) and 
the three configurations we consider being discontinuous, 
one could fear that a semi-classical analysis would fail. 
We see the relevance of the concept of Hawking tempera- 
ture as a confirmation that the configurations considered 
here are typical for observing the Hawking effect. In the 
next section we will draw the same conclusion from a 
study of the density correlations. 

From ([F7|) and Fig. [5] one gets an order of magnitude 



k^T^/{mcl) - 0.1, i.e., typically Th ~ 10 uK. Since the 
temperature in typical experiments is rather of the order 
of the chemical potential toc^ (i.e., around 100 nK), the 
Hawking radiation will be lost in the thermal noise and 
very difficult to identify. This is the reason why density 
correlations have been proposed in Ref. @ as a tool for 
identifying the Hawking effect. We thus consider two- 
body correlations in the next section. 



CORRELATIONS 



The connected two-body density matrix is defined by 



= {Xi,t)^\x2,t)^{xi,t)^{x2,t)) 

- {^\xi,t)^xi,t)){^Hx2,t)Mx2,t)). (58) 
g*-^^ is time independent because we work in a stationary 



configuration. In ([58)) the average is taken either on the 
groimd state or over a statistical ensemble, g'-^-' is directly 
related to the density correlations in the system, this can 
be seen by rewriting Eq. (|58p under the form 



g^'^\xi,X2) = {n{xi,t)n{x2,t)) ~ {n{xi)){n{x2)) 



(59) 



5{xi - X2){n{xi)). 



The last term in the right-hand side (r.h.s.) of (|59p is 
the Poissonian fluctuation term originating from the dis- 
creteness of the particles [1^. Written under this form, 
is 

sometimes denoted as the cluster function. 
For a system at thermal equilibrium in the grand 
canonical ensemble, one has 



n{x) = {n{x)) 
1 



= — Tr h{x) cxp 



fc„T 



(60) 



where Z = Tr{exp[-(i? - ^iN)/{k^T)]} is the partition 
function. Deriving expression ([50)) with respect to /i, one 
gets 



dn{x) 



{n{x)N) - {n{x)){N). 



(61) 



Since N = dx' n{x'), property (ICTI) and expression ([5^ 
yield the following sum rule: 



dx' g(^) {x, x') = -n{x) + k^T 



dn{x) 



(62) 



For a homogeneous system, this sum rule is a standard 
thermodynamic result (29j which can be shown to be 
equivalent to the compressibility sum rule (whose de- 
finition is given for instance in Ref. (soj). Formula 
(1621) is a generalization to inhomogeneous systems; it has 
been used in Ref. j31i] for witnessing quasi-condensation 
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through a study of density fluctuations and in Ref. (32| 
to propose an universal thermometry for quantum simu- 
lations. 

In the remainder of this seetion we concentrate on the 
T = case, postponing the discussion of finite tempera- 
ture to a future publication. Fulfillment of the sum rule 
([62|) is a strong test of the validity of the Bogoliubov ap- 
proach used in the present work. We give in Sec. IV Bl 
the leading order contributions to g^"^^ and explain in 
Appendix [C] how we use them in order to check that the 
T = version of the sum rule ()62p is indeed verified for 
— >■ cx) (i.e., far from the horizon). 

From the Bogoliubov expansion ([2]), one gets at leading 
order 

+ 'P{xi)'S>ix2){i^Hxi,t)i^\x2,t)) 

+ C.C.. (63) 



For i = 1 or 2, we write $(xi) — y/niexp(ikiXi)(l)i{xi), 
where m = n„ (n^), h = ku ikd) and 0i = 4>u {(pd) if Xi < 
{xi > 0). We recall that (j)d is defined in Eq. ^ and 0„ 
is either equal to unity (flat profile configuration) or given 
by Eq. (|TT|) (delta peak and waterfall configurations). 
Based on the decomposition (P5|) . one can show (see Ref. 
0) that Eq. §^ yields 



g'^^\xi,X2) = \/nin2 



2^ 



7(0:1, X2,W), 



(64) 



where 7(^1 , X2, uj) [and accordingly g^^\xi, X2)] is conve- 
niently separated in a zero temperature term 70 and a 
remaining term jj,: 



(65) 



-f{xi,X2,uj) = 7o(a;i,X2,a;) -I- 7r(a;i, 2:2, 1:^). 

f/Q and 7o are the contributions evaluated from ((63|) and 
(I45p which remain finite even in the T = case where 
HlIuj) = (ai(w)ai,(w)) = 0. One has 

'^o{xi,X2,i^) ^ ^ wl{xi)f^{x2) 

+ Q{n-uj)u*j^^{xi)fo2{x2)+c.c., (66) 

with 

UL{Xi) = (l)*{Xi)uLiXi 

and 

fi(a;j) = Mt(a;i) + Wt(a;i). 
The other contribution to (p5|) is 

7T(xi,a;2,a;) = f2(xi)fi(a;2)ni,(w) -|- c.c, (69) 



Wi,(a;i) = 0i(a;i)Wi(a;i), (67) 

(68) 



= E 



where it should be understood that the D2 contribution 
is only present for uj < fl. 



We often display below the results not for g^^^ but for 
the dimensionless quantity G^^' defined as 



G(2)( 



Xl,X2) 



g^^'>{xi,X2) 
n{xi)n{x2) 



(70) 



Also, we will compute the correlations when xi and X2 are 
both far from the horizon. In this case, 4ii{xi) = exp(i/3i) 
and G(2)(xi,a;2) = g'^'^\xi,X2) / {nin2). 



A. No black hole 

Before embarking in a general determination of g^"^^ for 
a black hole configuration, we recall here the result for a 
uniform fluid (density n„) moving at constant subsonic 
velocity. One gets from the no black hole version of ((66)) 

M 



jo{x,x',ll>) 



W£{x)ri{x') + c.c. 



E 



2TB7T 



£G{ti|in,u|out} 



(71) 



c.c. (72) 



This yields 



Gq {x, x') — 



where 



F{z) 



dt 



sin(2 1 z) 



(73) 



(74) 



(l-Ht2)3/2- 

This is the expected correlation in a quasi ID con- 
densate (see, e.g., Ref. [sH and references therein): 
nu^uG\) {x,x') is a universal function of z = (a; — x')/£^u- 
In particular, n„^„G[,^^(a;, x) = F{0) = -2/7r 36]. We 
evaluated the fluctuations around the uniform profile nu- 
merically by means of the truncated Wigner method for 
the Bose field [U HI] already used in [6| for studying the 
same observable in a black hole configuration. In Fig. 
[7] we compare the analytical form ([75]) with the results 
of the numerical computation along the cut displayed in 
the inset. The excellent agreement is a good test of the 
accuracy of the numerical method used in Ref. @ . 
We finally note here that dz F{z) = —1. This yields 



dx' g^^\x,x') = 



(75) 



which is a mere verification of the T = version of the 
sum rule (|62p in this simple uniform setting. 

The main correlation signal in the black hole configu- 
rations to be studied soon is similar to the short range 
anti-bunching displayed in Fig. [71 However, we will see 
that (i) its precise shape is affected in presence of an 
acoustic horizon and moreover, (ii) new long range cor- 
relations appear, which can be interpreted as emission of 
correlated phonons @. 
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FIG. 7. (Color online) C^n^G^^^ (x, x') computed analytically 
from Eq. (|73p (black solid line) compared with the numerics 
from the truncated Wigner method (red dashed line). The 
small discrepancy near a; — a;' = is due to numerical uncer- 
tainty and to the plotting procedure which introduces a small 
amount of smoothing of the raw numerical data. The inset 
represents a color plot of the numerical results. The white 
dashed line is the line x' = 50^i, — x along which we compare 
numerical and analytical results in the main plot. 



B. General formulae in presence of a black hole 



1. Case X and x' — >■ — oo 

We first consider the case where x and x' are both deep 
in the subsonic region, i.e., outside the black hole and far 
from the acoustic horizon. From Eq. (|66p . one gets in 
this case 



-fo{x,x',uj) = W*|in^«P'>e 



i9u|out(a:'-a;) 



+ W«|out'^n|oute- 

+ e(r! - c^)|5„,d2^7^„|outpe•«"l-'(^'"^) 

+ C.C.. (77) 

The contribution of the Su,d2 term disappears when cj > 
il and this is the reason for the Heaviside factor 9(ri — w) 
in ([77|) . If it were not for the Su,d2 term, ([T7)l would 
be exactly equal to ([7T|) . one would recover the same 
correlation as ([75)1 obtained in absence of black hole and 
the contribution of (|77p alone would be enough to verify 
the sum rule (|75p . Now the Su.d2 term is not zero, and 
this means that the correlations in the vicinity of the 
diagonal x = x' are modified by the existence of the 
black hole. This is similar to the results obtained by 
Kravtsov and coworkers for non standard ensembles of 
random matrices Indeed, for fixed x, ([77]) alone 

is not able to fulfill the sum rule. The addition of the non 
local correlations (1751) induced by the Hawking emission 
will be necessary to this end, as advocated in Ref . ^] . 



We now turn to the study of zero temperature density 
fiuctuations around a sonic horizon. For simplicity, we 
only consider the case where x and x' are far from the 
horizon: this makes it possible (i) to drop the evanescent 
contributions to (j66P and (ii) to avoid treating the posi- 
tion dependence of the background density in the delta 
peak and waterfall configurations. Part of these results 
were already obtained in Ref. @ (the ones valid when x 
and x' are far one from the ot her) and here we generali- 
ze and correct some misprints [39| . We only display the 
most important contributions to 70 which are both the 
larger ones and the ones useful for fulfilling the T = ver- 
sion of the sum rule (I62P when x is far from the horizon 
[this version is written explicitly in Eq. (|8ip]. We note 
here that a similar approach has previously been used in 
Ref. [S^l for studying phase fluctuations in a similar set- 
ting, with a description of the scattering less elaborate 
than the one presented in Sec. IIIIEI 

Finally, we introduce the notations Ue and We defined 
by 

\m{x))-'' \m{x))' ^^^^ 

where ut{x) and wt,{x) are defined in Eqs. ([S7)) . We also 
introduce TZi =Ue + VV^. 



2. Cases (x — > —00 and x' — > +00^ or (x +00 and 
x —00 ) 

In the case where x is deep in the upstream region and 
x' deep in the downstream one, we get 

7o(x,2:',w) = e(17-a;)S',*,i2'5'di,d2X 
+ e(n-u;)S;rf2^d2,d2X 

'^ti|out '^ci2|out"= 

-l-cc. (78) 

If instead x is deep in the downstream supersonic region 
and x' deep in the upstream subsonic region, it suffices 
to exchange the roles of x and x' in the above formula. 

3. Case x and x' +00 

This is the case where x and x' are both deep in the 

downstream region (i.e., deep inside the black hole). The 

(2) 

leading order contribution to can be separated in a 
diagonal part which depends only on x — x' and a non 
diagonal part. The diagonal part reads 

jtHx,x\uj) = W:iM„^di|i„e''-i'"(-'--) 
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d2|out '^(i2|out 
+ e(f}-^)(|S'rf2,u|' + |^d2,dinX 

|7^d2|out|'e''-|-(^'-^) 
+ C.C.. (79) 

In absence of black hole, the terms involving coefficients 
of the S'-matrix disappear in ([7^ and this gives, after 
integration over w £ M"*" , the usual quasi-condensate cor- 
relation signal: '{x,x') = ndF[{x - x')/£,d\/£.d- 

The non diagonal part is only present if a horizon exists 
and only contributes for w < il; it reads 



l0°'^~'^"^^(x,x',Uj) = Q{n- Uj)S2i^d2Sd2,d2^ 
'^£;i|out'^d2|out^' 

+ {x< — ^x')+c.c.. (80) 



C. Results for the three configurations 

Formulae ^ and §UI allow us to deter- 

(2) _____ 

mine {x,x') through Eq. (|64l) . We performed the 
corresponding integration over uj G M"*" numerically. The 
results are shown in Fig. [5] for the delta peak configu- 
ration and in Fig. [9] for the waterfall configuration. In 
each of these figures we only display the correlations for 
I a; I and \x'\ larger than a few healing lengths, because we 
use formulae which are exact only in the limit |a;| and 
\x'\ -5> oo. 

In each plot, we also display the lines where the heuris- 
tic interpretation of the Hawking effect presented in the 
introduction (see also Ref. ^) leads to locate the more 
pronounced correlation signal: if correlated Hawking 
phonons are emitted along the wjout, dljout and c?2|out 
channels, at time t after their emission, these phonons 
are respectively located at a;„|out — (Ki — Cii)t < 0, 
^i|out = {Vd + Cd)t > and a;d2|out = {Vd - Cd)t > 
[40]. This induces a correlation signal along lines of 
slopes: {Vu — c„)/(Vd -I- Cd) (resulting from u — dl cor- 
relations), {Vu — Cu)/{Vd — Cd) {u — d2 correlations) and 
{Vd — Cd)/{Vd + Cd) {dl — d2 correlations). Of course the 
lines with inverse slopes are also present (they correspond 
to the exchange x O x'). Indeed, the main features of 

f 2) 

the computed perfectly match this interpretation of 
the Hawking effect. 

These results are very similar to the ones obtained nu- 
merically for the flat profile configuration (already dis- 
played in Refs. @,|1]). This legitimizes the use of density 
correlations as a tool for identifying Hawking radiation 
also in the realistic delta peak and waterfall configura- 
tions. 
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FIG. 8. (Color online) 2D plot of the numerical result for the 
quantity ^uriuGg^' {x, x') in the case of a delta peak configura- 
tion with m„ = 0.5. The shaded area near the axis corresponds 
to the zone \x\ or \x'\ < lO^u- Gg^' is only displayed for 
and \x'\ > lO^u, i-e., in the asymptotic region where expres- 
sions JTtIi, ((ZSl), ([79} and ((80]) are valid. The dashed straight 
lines correspond to the correlation lines where a heuristic in- 
terpretation of the Hawking signal leads to expect the largest 
signal (see the text). 
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FIG. 9. (Color online) Same as Fig. [8] for a waterfall config- 
uration with fill, = 0.5. 



For each plot the dominant signal is the anti-bunching 
along the diagonal {x — x'). This corresponds to the 
typical local density correlation in a quasi-condensate (cf. 
Fig-IZl)- However, this signal is modified compared to the 
one observed in a uniform system (see, e.g., the discussion 
in Sec. IVB 1|) . This modification is connected to non 
local features which are necessary to verify the sum rule 



dx' gj^\x, x') = 



-riu when x 
-Ud when x 



-oo, 
-oo, 



(81) 



which is the T — version of the sum rule (p^ valid 
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when X is far from the horizon. We checked this sum rule 
analytically in Appendix [Qon the basis of our Bogoliubov 
description of the quantum fluctuations and of the results 
of Sees. IVBTI IV¥2l and IVB3l 

By comparing Figs. [8] and |9] with the inset of Fig. 
[7l one can reverse the argument of Ref. [12] and ar- 
gue that (i) in presence of an acoustic horizon, the main 

new features of are the non local density correlations 
which are simply understood as resulting from the emis- 
sion of correlated phonons and (ii) because of the sum 
rule (I5T|) . these long range correlations have to be associ- 

(2) 

ated to modifications of the short range behavior . 
However, these short range modifications are not of great 
experimental relevance because they would be efficiently 
blurred by finite temperature effects. The non local as- 
pects instead are good signatures of the Hawking effect 
because they are easily distinguished from the thermal 
noise (being even reinforced at finite T as demonstrated 
in Ref. [|). 

VI. CONCLUSION AND DISCUSSION 

In the present work we have introduced new and realis- 
tic acoustic analogs of black holes and have analyzed the 
associated Hawking radiation. We restricted ourself to 
simple configurations (flat profile, delta peak, waterfall) 
but our approach is easily applied to more complicated 
cases. For instance, in a wave guide with a constriction, 
one is in a mixed situation where there is an external 
potential step (as in the waterfall configuration) whereas 
the non linear parameter is position-dependent (as in the 
flat profile configuration) [13|. This case is interesting 
since it may be possible to realize it experimentally, but 
the analytical treatment is straightforward and we do 
not consider it here because this would bring no new in- 
sight on our theoretical method. One could also consider 
smooth potentials, and the eigen-modes (defined in Sec. 
imp should then be determined numerically, but the the- 
oretical framework presented here remains of course valid 
in this situation. We note also that the present approach 
can be adapted to treat the creation of a black hole hori- 
zon in a Fermi gas as suggested in Ref. [4l| . 

The description of the system in terms of a S'-matrix 
allows for a simple description of the radiation spectrum 
and a clear identification of the characteristics of the 
system. In particular, we showed that Hawking radia- 
tion is absent for a "no black hole configuration" corre- 
sponding to a flow connecting two subsonic asymptotic 
regions. The spectrum of Hawking radiation has been 
computed in the case of a black hole connecting an up- 
stream subsonic region to a downstream supersonic one, 
and the concept of Hawking temperature has been dis- 
cussed quantitatively. 

The main focus has been put on non local density cor- 
relations. We verified that their interpretation in terms 
of emission of correlated phonons previously introduced 
in a model configuration [6i, i9j] also holds in realistic set- 



tings. By studying a sum rule verified by the two-body 
density matrix, we showed that the Bogoliubov descrip- 
tion of the quantum fluctuations around the stationary 
ground state of the system also provides an accurate de- 
scription of the short range density correlations. 

In this work we introduced new acoustic black hole 
configurations motivated by their possible experimental 
realization and proposed - following Ref. @ - non local 
density correlations as practical signatures of Hawking 
radiation. It is thus important to discuss if the proposed 
signal is large enough for being detected experimentally. 
From Figs. |S]and[ni one sees that the prominent Hawking 
signal corresponds io u — d2 correlations. For each figure 
this corresponds to a line of negative correlation where 
the largest value of ^uTIuGq is between —0.01 and —0.02. 
In present-day experiments it is possible to measure den- 
sity fluctuations around a mean density Uu of order 10 
/im^^ in a setting where the healing length ^„ is around 
a few ^m (see, e.g., Ref. [3l|). It is thus realistic to hope 
to reach a configuration where |Gg^^|max — 5 x 10"-^, 
and for detecting a signal of this intensity a precision of 
around 10~'* is required on the determination of G^^'. 
Noticing that G*-^-' is the quadratic relative density fluc- 
tuations, detecting this signal would correspond to mea- 
suring density fluctuations with a precision of order of 1 
%, which seems within reach of present-day experimental 
techniques. 
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Appendix A: Low energy behavior of the scattering 
matrix 

In this appendix we display the analytical results for 
the low-w behavior of combinations of the elements of 
the S'-matrix relevant for computations of the Hawking 
temperature and for the fulflllment of the sum rule (j8ip . 
We only give results for the flat proflle and waterfall con- 
figurations because those concerning the delta peak con- 
figuration are too long. Indeed, for the delta peak con- 
figuration, a numerical determination of the coefficients 
of the 5-matrix (which simply amounts to invert a 4 x 4 
matrix) is more convenient than the analytical approach. 
We checked that both agreed to an extremely good ac- 
curacy. 
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1. Flat profile configuration 



In this case, the scattermg coefficients depend of the 
two Mach numbers m„ and md- The coefficients / and h 
defined in Eqs. (H^ verify the following relations 



|/u,(i2p — 2- 



"td "t^ — 1 + ntu 



(Al) 



\Jdl,d2\ — 77 



\fd2,d2\ — t: — 



fu,d2f< 



dl,d2 



fu,d2fd2,d2 



fdl.d2fd2, 



md 


+ mu 


1 + 


mu 


md 


+ mu 


1 - 


mu 


md 


— mu 


1 + 


mu 




mu 


1 - 


mu 


md 


+ mu 


14 


' mu 


3 








2 


mu 


1 


— m 



{mj-l)K (A2) 



(mj-l)^, (A3) 



{mj~ir^, (A4) 



md J md — m 



u 1 + ntw 



imj-l)-2, (A5) 



d2 



md 



1 



Re(/* W2^u,d2) 



(l + -«)2- 

2. Waterfall configuration 



(mi-1)^ (A6) 



(A7) 



fu,d2}d2.d2 



^i{l + mu)^{l-muY^{l + inly. 



(1 



m-a + ^2)2 



(A12) 



. _ 1 {l + mu)^{l-mu)i{l + ml)i 
Jdl,d2Jd2,d2 - 2 ~, ; 2V2 ' l^-'^'^J 



R-e(/*.d2^«, 



d2) 



{l + mu + mlY 
«m(2 + '«^)(1 + 2m2) 



(A14) 



(1 + ;„„)2(1 + ;„„ + ;„2)2 

Appendix B: Computation of the energy current 

In this appendix we briefly indicate how formula (j53p 
is obtained from Eqs. (HH]) and ((5(I| . For a point x 
deep in the subsonic region, using the ry-unitarity (j43p 
of the S'-matrix, one can write the relevant contributions 
to Jo(a;,i^) in (|50p under the form 

Jo(a;,w) = (g„|i„ - A:„)|W„|inP + (g„|out - ^«)|W„|outr 



d2\ 



(9ti|out + fciOl^uloutl^ 
f (gnlout - ^«)|>Vn| out I 



(Bl) 



A simple but lengthly computation shows that the con- 
tributions to Ho of the two first terms of the r.h.s. of (jBip 
cancel after integration over w. This is very satisfactory 
because this shows that there is no Hawking radiation 
when Su.d2 = 0, i.e., in absence of black hole. 

Using Eq. (j25p the remaining can be written as 
\Su,d2\'^Ju\out/^uCu, and since h/{mCuS,u) = 1 this directly 
yields Eq. (|53p . The minus sign in this formula comes 
from the fact that J„|out = ~1 corresponds to the 
direction of propagation of the energy in the u|out mode. 



In the case of the waterfall configuration, the elements 
of the 5-matrix depend on a unique parameter; we chose 
to express them as functions of the Mach number mu ■ 



|/M,d2p — 2 



'"m(1 - m„)2 (1 + ml)2 
(l + m„)^(l + m„ + m2)2 



\fdl,d2 



2_1 {l-Mu)--[l + ml)i 



2{l + mu)"-{l + mu + mlY' 



I/. 



d2,d2\ 



1 {^~mt)i 
2{l + mu + mlY' 



/u,d2/rfl,rf2 



m^jl - m„)2(l + ^2)2 
(l + m„)5(l + m„ + m2)2 



(A8) 



(A9) 



(AlO) 



(All) 



Appendix C: Verification of the sum rule (|62|) 

In this appendix we check that the Bogoliubov ap- 
proach based on expansion (j45p indeed makes it possi- 
ble to verify Eq. ([5T|) which is the T = version of the 
sum rule (|62p when x is far from the sonic horizon, either 
upstream or downstream. 

We start here by a technical remark. For fixed x, the 
contributions to ^q{x^x' ,uj) displayed in Sees. IV B 1[ 
IV B 21 and IVB 31 are only noticeable for x' ^ — ^„ or 
x' ^ £,d- As a result, the analytical forms displayed 
in these sections can be extended for all G M without 
introducing noticeable errors in the computation of the 
integral of "fo{x, x' ,uj) over x' . Then, the integration of 
7o(x, x' ,Lo) over x' just amounts to evaluating the follow- 
ing integral: 



(CI) 
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J if ^ e {dl|in,d2|in}, 

I iVgiQei^Wi^^) if ^ e {u|in,M|out,dl|out,(i2|out}. 

In (|C1|) we used the fact that qdi\in and 9^2 1 in never can- 
cel, whereas the other g^'s do for w = (see Fig. [3]). 
Using this prescription, for x large and negative, one gets 
from Eq. 



dx' Tiu I ^7o(x,x',w) 

ZTT 



I fu 



42\ 



2 1 



(C2) 



and from Eq. ([751) 

dx' y/nuUd 







_ n« / Cu rid -p^^ 
2 V Cd n„ 

Altogether this yields 



doj 
2^ 

/* 
u,d2 



1 



7o(x,x',a;) 

{fdl,d2 + fd2,d2) 



(C3) 



lim / dx'gQ{x,x') 



2 V V 1 



J u. 



(12 



(C4) 



where 



— fu,d2\ —— + fdl.d2 + fd2.d2- (C5) 
V C„ rid 

Similarly one gets [from Eqs. dTS]), ^ and dHHl)] 



lim 




x—>--\-oo 










.Cd 



-rid 



Idl,d2 



I d2A2 



md 



1 



md 



T 



(C6) 



Using the analytical expressions for the combinations of 
coefficients /„,rf2, Jd\,d2 and /d2,d2 displayed in Appendix 
1X1 one can easily verify that the second terms of the 
r.h.s. of Eqs. (jC4l) and (IC6p cancel. This is due to the 
fact that J- is identically null in the flat profile and water- 
fall configurations (this is more tedious to check but we 
confirmed it analytically). The same holds for the delta 
peak configuration. This shows that the sum rule ((5T|) is 
fulfilled in these three cases. This is a strong confirma- 
tion of both the validity of the Bogoliubov approach and 
of the exactness of our analytical results. 
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